// P21Diffractor.cpp: implementation of the P21Diffractor class.
//
//////////////////////////////////////////////////////////////////////


#include "P21Diffractor.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

///////////////////////////////////
// b unique
///////////////////////////////////
P21Diffractor::P21Diffractor(char *DiffractorFile, 
						   Position *W1, Position *W2, Position *W3, char **atomLabel, int nAtom)
: Diffractor(DiffractorFile, W1, W2, W3, atomLabel, nAtom)
{
newA(Case, NIntensity);
newA(a1, NIntensity);
newA(a2, NIntensity);
newA(a3, NIntensity);
newA(b1, NIntensity);
newA(b2, NIntensity);
newA(b3, NIntensity);
newA(AP, NAtom);

for(int n=1;n<=NIntensity;n++)
	Case[n] = ( K[n]%2 == 0 );
}


P21Diffractor::~P21Diffractor()
{
deleteA(Case);
deleteA(a1);
deleteA(a2);
deleteA(a3);
deleteA(b1);
deleteA(b2);
deleteA(b3);
deleteA(AP);
}


//////////////////////////////////////////////////////////////////////
// functions
//////////////////////////////////////////////////////////////////////

void P21Diffractor::TransformDerivatives(double   *IDerivative, 
										 Position *XDerivative)
{
int    n, m, h, k, l;
double c1,c2,s1,s2;
double x1,x2;
double F1,F2,F3,*FF;

for(n=1;n<=NIntensity;n++)
	{
	x1  = IDerivative[n]*ReAmplitude[n];
	x2  = IDerivative[n]*ImAmplitude[n];

	a1[n] = x1*DH[n];
	a2[n] = x1*DK[n];
	a3[n] = x1*DL[n];

	b1[n] = x2*DH[n];
	b2[n] = x2*DK[n];
	b3[n] = x2*DL[n];
	}


for(m=1;m<=NAtom;m++)
	{
	SetUpTrig(AP[m]);

	FF = FiniteFactor[m];

	F1 = 0.0;
	F2 = 0.0;
	F3 = 0.0;

	for(n=1;n<=NIntensity;n++)
		{
		h = H[n];
		k = K[n];
		l = L[n];

		s1 = (SinX[h]*CosZ[l] + CosX[h]*SinZ[l]);
		c1 = (CosX[h]*CosZ[l] - SinX[h]*SinZ[l]);
		c2 = FF[n]*CosY[k];
		s2 = FF[n]*SinY[k];

		if(Case[n])
			{
			F1 -= s1*(c2*a1[n] + s2*b1[n]);
			F2 += c1*(c2*b2[n] - s2*a2[n]);
			F3 -= s1*(c2*a3[n] + s2*b3[n]);
			}
		else
			{
			F1 += c1*(c2*b1[n] - s2*a1[n]);
			F2 -= s1*(c2*a2[n] + s2*b2[n]);
			F3 += c1*(c2*b3[n] - s2*a3[n]);
			}
		}

	FDeriv[m][1] = F1;
	FDeriv[m][2] = F2;
	FDeriv[m][3] = F3;
	}
			

MakeXDerivative(XDerivative);	
}



void P21Diffractor::FourierTransform(Position *AtomPosition)
{
int    n, m, h, k, l;
double c, s, *FF;

for(m=1;m<=NAtom;m++)
	AP[m] = AtomPosition[m];

for(n=1;n<=NIntensity;n++)
	{
	ReAmplitude[n] = 0.0;	
	ImAmplitude[n] = 0.0;
	}


for(m=1;m<=NAtom;m++)
	{
	SetUpTrig(AtomPosition[m]);

	FF = FiniteFactor[m];

	for(n=1;n<=NIntensity;n++)
		{
		h = H[n];
		k = K[n];
		l = L[n];

		if(Case[n])
			{
			c = FF[n]*(CosX[h]*CosZ[l] - SinX[h]*SinZ[l]);

			ReAmplitude[n] += c*CosY[k];
			ImAmplitude[n] += c*SinY[k];
			}
		else
			{
			s = FF[n]*(SinX[h]*CosZ[l] + CosX[h]*SinZ[l]);

			ReAmplitude[n] -= s*SinY[k];
			ImAmplitude[n] += s*CosY[k];
			}
		}
	}
}